
clear all
clc

load('Match_EndShare_minmax_small.mat')

b = str2double(getenv('SLURM_ARRAY_TASK_ID'));
%b = 1;

hat = 0; hat_ind = 0;
bonus_hat = 1;
v       = 2;

%% Model parameters

datsize = 100;
start   = 2;
S       = 200;

% Data samples
if b == 1
    sample = 13;
    keep = find(data(:,sample) == 1);
    data = data(keep,:);
    keep = find(data0(:,sample) == 1);
    data0 = data0(keep,:);
    clear keep
end

% Imputted bonuses
if bonus_hat == 1
    bonus_ind = 1;
else
    bonus_ind = 0;
end

%% Define data matricies and parameter counts

N = size(data,1);
propID = data(:,1);
firmID = data(:,5);

% Data needs to be sorted by market
market = data0(:,2);
marketindex = unique(market);
marketindex = [(1:size(marketindex,1))' marketindex];

M = max(marketindex(:,1));
MarketSizeVec = zeros(M,1);
for m = 1:M
    row = find(market == marketindex(m,2));
    MarketSizeVec(m,1) = size(row,1);
end
J = unique(data(:,5)); J = size(J,1);

FirmMatch = data(:,5);
FirmMatch0 = data0(:,5);

HOLD = [data0(:,5) data0(:,8)];
FirmMatchConstraint = [(1:J)' zeros(J,1)];
for j = 1:J
   row = find(HOLD(:,1) == j);
   FirmMatchConstraint(j,2) = size(row,1);
end
clear HOLD
FirmMatchConstraint = FirmMatchConstraint(:,2);

%% Share data

% Stacked share of leases signed by a firm in a wellpad
share0 = data(:,7)./data(:,43);
share0(isnan(share0)) = 0;

% Firm share observed in the data to match in the algorithm
FirmShare0 = [data(:,2) data(:,5) data(:,7) data(:,43)];
FirmShare0 = unique(FirmShare0,'rows');
FirmShare0 = FirmShare0(:,3)./FirmShare0(:,4);
FirmShare0(isnan(FirmShare0)) = 0;

%% Data vectors

pipeDist        = data(:,49); pipeDist0        = data0(:,49);
wellDist        = data(:,50); wellDist0        = data0(:,50);
landSize        = data(:,48); landSize0        = data0(:,48);
prodwell        = data(:,27);  prodwell0       = data0(:,27);
enforce         = data(:,28);  enforce0        = data0(:,28);
app_val         = data(:,17);  app_val0        = data0(:,17);
parcelToWP      = data(:,57);  parcelToWP0     = data0(:,57);
parcelToApi     = data(:,60);  parcelToApi0    = data0(:,60);
parcelInt       = data(:,16);  parcelInt0      = data0(:,16);
landman         = data(:,33);  landman0        = data0(:,33);
largeOp         = data(:,21);  largeOp0        = data0(:,21);

% Lease quality
Clause  = data(:,67 + hat_ind); Clause0 = data0(:,67 + hat_ind); 

% Well profit
expFutureParcel = data(:,53 + bonus_ind);  expFutureParcel0  = data0(:,53 + bonus_ind);
expFutureFirm   = data(:,51 + bonus_ind);  expFutureFirm0  = data0(:,51 + bonus_ind);

firmcount = data(:,11); firmcount0 = data0(:,11);

ParcelIndex     = [expFutureFirm parcelToWP landSize app_val pipeDist wellDist parcelInt Clause];
Operator        = [landSize.*largeOp wellDist.*largeOp parcelInt.*largeOp];
Landman         = [landSize.*landman wellDist.*landman parcelInt.*landman];

ParcelIndex     = [ParcelIndex Landman Operator];

ParcelIndex0    = [expFutureFirm0 parcelToWP0 landSize0 app_val0 pipeDist0 wellDist0 parcelInt0 Clause0];
Operator0       = [landSize0.*largeOp0 wellDist0.*largeOp0 parcelInt0.*largeOp0];
Landman0        = [landSize0.*landman0 wellDist0.*landman0 parcelInt0.*landman0];

ParcelIndex0    = [ParcelIndex0 Landman0 Operator0];

FirmMatrix      = [expFutureParcel  prodwell  enforce  wellDist.*largeOp   landSize.*largeOp  Clause];
FirmMatrix0     = [expFutureParcel0 prodwell0 enforce0 wellDist0.*largeOp0 landSize0.*largeOp0 Clause0];

NN = size(ParcelIndex,2);
JJ = size(FirmMatrix,2);

%% Halton draws (use a different set of S draws for each MC)
B = v;
StateVec = randi(1e9,[B,2,S]);

FirmError_Sim = zeros(N,S);
ParcelError_Sim = zeros(N,S);

for s = 1:S

const.NObs = N; 

rng(StateVec(v,1,s))
skip = randi(1e3);
leap = randi(1e2);

p = haltonset(1,'Skip',skip,'Leap',leap);
p = scramble(p,'RR2');
draw = net(p,const.NObs);
draw = norminv(draw);

% FirmError_Sim(:,s) = draw(firmID,:);
FirmError_Sim(:,s) = draw;

clear draw p

const.NObs = N; 

rng(StateVec(v,2,s))
skip = randi(1e3);
leap = randi(1e2);

p = haltonset(1,'Skip',skip,'Leap',leap);
p = scramble(p,'RR2');
draw = net(p,const.NObs);
draw = norminv(draw);

ParcelError_Sim(:,s) = draw;

clear draw p

end

clear skip leap

display(['FirmError: ' num2str(sum(sum(FirmError_Sim))) ])
display(['ParcelError: ' num2str(sum(sum(ParcelError_Sim))) ])

%% Market Index

MarketCount = zeros(M,5);
MarketIndex = zeros(1,M);
begin = 1;

for mm = 1:M
    [row] = find(data(:,2) == marketindex(mm,2));
    [row0] = find(data0(:,2) == marketindex(mm,2));

    firm = data(row,5); firm = unique(firm);
    MarketIndex(1:size(row,1),mm) = row;

    MarketCount(mm,1) = size(row,1);  % N for each market
    MarketCount(mm,2) = begin;        % Beginning index for each market
    MarketCount(mm,3) = MarketCount(mm,2) + size(row,1) - 1; % Ending
    MarketCount(mm,4) = size(firm,1); % J for each market
    MarketCount(mm,5) = MarketCount(mm,1)/size(row0,1);
     
    begin = begin + size(row,1);
end

% 1: count of land in market; 2: beginning index; 3: ending index
MarketCount0 = zeros(M,5); begin = 1; beginST = 1;
for mm = 1:M
    [row] = find(data0(:,2) == marketindex(mm,2));
    MarketCount0(mm,1) = size(row,1);  % N for each market
    MarketCount0(mm,2) = begin;        % Beginning index for each market
    MarketCount0(mm,3) = MarketCount0(mm,2) + size(row,1) - 1; % Ending

    begin = begin + size(row,1);
    
    [row] = find(data(:,2) == marketindex(mm,2));
    
    MarketCount0(mm,4) = beginST;        % Beginning index for each market
    MarketCount0(mm,5) = MarketCount0(mm,4) + size(row,1) - 1; % Ending

    beginST = beginST + size(row,1);
end

%% Data for the observed match

MarketIndex0 = zeros(1,M);
index = 1;

for mm = 1:M
    [row] = find(data0(:,2) == marketindex(mm,2));
    MarketIndex0(1:size(row,1),mm) = [index:index + size(row,1) - 1]';
    index = index + size(row,1);
end

clear begin index mm row

%% Prepare Estimation
% gamma = ones(3,NN-1)*0.1; delta = ones(3,JJ-1)*0.1;
gamma = rand([3,NN - 1]); delta = rand([3,JJ - 1]); 
alpha = rand([3,1]);

theta = [gamma(start,:) delta(start,:) alpha(start,:)];

%% Moments for the data

match.ParcelError_Sim   = ParcelError_Sim;
match.FirmError_Sim     = FirmError_Sim;

match.MarketIndex       = MarketIndex;

match.MarketCount       = MarketCount;
match.FirmMatch         = FirmMatch;
match.ParcelIndex       = ParcelIndex;
match.FirmMatrix        = FirmMatrix;

match.MarketIndex0      = MarketIndex0;
match.FirmMatch0        = FirmMatch0;
match.ParcelIndex0      = ParcelIndex0;
match.FirmMatrix0       = FirmMatrix0;

match.NN                = NN;
match.JJ                = JJ;
match.N                 = N;
match.J                 = J;
match.M                 = M;
match.S                 = S;

match.FirmShare0        = FirmShare0;             % Share by market (J*M x 1)
match.share0            = share0;                 % Stacked share (J*N x 1)

match.FirmMatchConstraint   = FirmMatchConstraint;
match.MarketCount0          = MarketCount0;

MomentSwitch = [1 1 1 1 1 ];
match.MomentSwitch = MomentSwitch;

[PMean0,FMean0,JntDist0,Cov0,Pvar0,Fvar0,PFcovar0,Pmean_firm0,Pvar_firm0,PFJntDist_firm0,PFcovar_firm0,GpMean0,GpVar0,GpMean_byJ0,GpVar_byJ0] = GroupMoments0_new(match.ParcelIndex0,match.FirmMatrix0,match.FirmMatch0,J,MarketIndex0,MarketCount0,MomentSwitch);
[GpPMeanShare0,GpFMeanShare0,GpShareMarket0,GpShareFirm0,GpFMean0,GpParcelCov0,GpFirmCov0] = GroupMoments0_share2(match.ParcelIndex0,match.FirmMatrix0,match.FirmMatch0,J,MarketIndex0,MarketCount0,FirmShare0);

temp                    = [FMean0 PMean0];
match.MMean_Moment      = temp(:);
match.MCov_Moment       = PFcovar0(:);
temp                    = [Fvar0 Pvar0];
match.MVar_Moment       = temp(:);

match.JntDist0          = JntDist0(:);
match.Cov0              = Cov0(:);

match.JMean_Moment      = Pmean_firm0(:);
match.JCov_Moment       = PFcovar_firm0(:);
match.JVar_Moment       = Pvar_firm0(:);
match.JDist_Moment      = PFJntDist_firm0(:);

% match.meanGroup_Moment  = GpMean0(:);
% match.varGroup_Moment   = GpVar0(:);

match.JMVar_Moment      = GpVar_byJ0(:);
match.JMMean_Moment     = GpMean_byJ0(:);

temp                    = [GpFMeanShare0 GpPMeanShare0];
match.GpMeanShare0      = temp(:);
match.GpShareMarket0    = GpShareMarket0;
match.GpShareFirm0      = GpShareFirm0;
match.ShareFMean0       = GpFMean0;
temp                    = [GpFirmCov0 GpParcelCov0];
match.GpCovShare0       = temp(:);

%% Estimate

tic
gmm_myopic = mindistGrpOfferGS_out_myopic2(theta,match)
toc

time = clock;

% Genetic algorithm
gaoptions = gaoptimset('MutationFcn',{@mutationgaussian, 4, 0.25}, ...
       'Generations',75,...
       'Display','iter',...
       'CrossoverFraction',0.75,...
       'EliteCount',1,...
       'PopulationSize',40);
[betaStar,fval,exitflag,output,population,scores] = ga(@(theta)mindistGrpOfferGS_out_myopic2(theta,match),size(theta,2),[],[],[],[],[],[],[],gaoptions);

% Pattern search
% patternoptions = psoptimset('MaxIter',50,'cache','on','Display','iter','ScaleMesh','off','SearchMethod',@positivebasisnp1);
% % ,'UseParallel','true'
% [betaStar,fval,exitflag,output] = patternsearch(@(theta)minimumdistance(theta,match),theta,[],[],[],[],[],[],[],patternoptions);

time = etime(clock,time)

betaStar

clear ParcelError_Sim FirmError_Sim

% Fminsearch
options = optimset('Display','iter','Diagnostics','off','TolX',1e-8,'TolFun',1e-8,'MaxIter',1e5,'MaxFunEvals',1e6,'LargeScale','off');
[theta_score,fval,exitflag] = fminsearch(@(theta)mindistGrpOfferGS_out_myopic2(betaStar,match),betaStar,options);

theta_score

% Local search with nelder meade

% tau = 100*ones(size(betaStar,2),1);
% ns = 500;		   
% nt = max(50,10*size(betaStar,2));
% eps = 1e-9;
% n = size(betaStar,2);
% 
% time = clock;
% [xmin,ynewlo,icount,numres,ifault] = nelmin ( @(betaStar)mindistGrpOfferGS_out_myopic(betaStar,match), n, betaStar,eps,tau,ns,nt);
% time = etime(clock,time)
% 
% xmin 

clear ParcelError_Sim FirmError_Sim
clear match

save(['DataResult_v' num2str(v) '_ival' num2str(b) '_myopic_out_' num2str(datsize) '.mat' ], '-v7.3')

quit

